Agenda
Introducción al modelado de regresión
Modelo de Regresión Lineal Simple
Regresión Lineal Simple en R
Regresión Lineal Múltiple
Regresión Lineal Múltiple en R
Evaluación de Supuestos
El Modelo Lineal Generalizado
La regresión de Poisson
La regresión Binomial Negativa
Tablas de regresión reproducibles con {gtsummary}
Conjunto de técnicas estadística para estimar la relación entre variables.
Los modelos de regresión multivariable modelan una sola variable dependiente en función de una o más variables independientes.
Según STRATOS podemos usar regresión para 3 propósitos diferentes:
Descripción*
Predicción
Explicación
Este curso se centrará solamente en algunas aplicaciones.
No abordaremos modelos de regresión para desarrollar modelos o reglas de predicción clínica.
Tampoco para métodos de inferencia causal robusta.
Agenda
Introducción al modelado de regresión
Modelo de Regresión Lineal Simple
Regresión Lineal Simple en R
Regresión Lineal Múltiple
Regresión Lineal Múltiple en R
Evaluación de Supuestos
El Modelo Lineal Generalizado
La regresión de Poisson
La regresión Binomial Negativa
Tablas de regresión reproducibles con {gtsummary}
Método estadístico que modela la
relaciónentre unavariable continua (dependiente)y otrasvariables (independientes).
\(Y\) es variable resultado (outcome), respuesta o dependiente.
\(X\) es una variable explicativa, predictora o regresora.
En la figura, a mayor valor de \(X\), mayor valor de \(Y\).
Podemos tratar de dibujar una línea recta que resuma la relación.
Existen infinitas rectas posibles que podríamos trazar: ¿Cuál elegir?
recta que pase por el valor más representativo del \(y_i\) en cada valor fijo de \(x_1\).
recta que conecte los promedios condicionados en \(x_1\)regresión lineal simple se puede expresar de la siguiente manera:\[y_i = \underbrace{\beta_0 + \beta_1x_{1i}}_{\text{componente} \\ \text{sistemático}} + \underbrace{\epsilon_i}_{\text{componente} \\ \text{aleatorio}}\]
El componente sistemático describe la media de \(y_i\) en cada valor fijo de \(x_{1i}\): Media condicionada de y en x.
El componente aleatorio describe la variación de los individuos alrededor de cada media condiciona de y:
observación \(i\) en la población, podemos relacionar el valor esperado (media) \(E[y_i]\) de \(y_i\) (también llamado \(\mu_i\)) con la variable explicativa \(x_{1i}\) mediante la siguiente ecuación lineal:\[E[Y | X_1 = x_{1i}] = E[y_i] = \mu_i = \beta_0 + \beta_1x_{1i}\]
variables aleatorias independientes e idénticamente distribuidas (i.i.d)parámetros desconocidos de una superpoblación infinita.\(x_1\) es fijo
\(\beta_0\) y \(\beta_1\)
coeficientes de regresión y son una medida de asociación.queremos estimar con los datos de la muestra!Advertencia
componente sistemático solo relaciona el promedio condicionado de \(y_i\) con las variables explicativas, NO con los valores individuales.\[\epsilon_i = y_i - \mu_i\]
El problema es que el término de error \(\epsilon_i\) no puede predecirse ni estimarse con los datos, se considera que es el componente no explicado por la variable independiente.
\[y_i = \beta_0 + \beta_1x_{1i} + \epsilon_i\]
\[y_i \sim N(\beta_0 + \beta_1x_{1i}, \sigma^2)\]
\[\epsilon_i \sim N(0, \sigma^2) \]
Usamos métodos numéricos:
Método de Mínimos Cuadrados Ordinarios (MCO)
Método de Máxima Verosimilitud (MV)
MCO y MV son equivalentes para el caso de la regresión lineal normal.
El estimador MCO es insesgado, no importa la distribución de \(y_i\) o \(\epsilon_i\).
El estimador MCO tiene mínima varianza si y solo si:
independencia de observaciones.normalidad.No es necesario que \(\epsilon_i\) o sigan una distribución normal para que los coeficientes de regresión \(\beta\) puedan estimarse de manera puntual.
Sin embargo, para estimar el valor p o los intervalos de confianza mediante inferencia clásica sí se necesita asumir una distribución conocida.
El modelo de regresión lineal normal asume normalidad de \(y_i\) y \(\epsilon_i\).
Sin embargo, el modelo es robusto a desviaciones leves/moderadas de la normalidad cuando se cumple el TLC (número de observaciones grande).
Otros enfoques para inferencia flexibilizan este supuesto:
Las variables categóricas no son continuas, en cambio son discretas y asumen solo unos cuantos valores.
¿Cómo estimar una medida de asociación cuando la variable explicativa es categórica?
Veamos el caso binario.
Si la variable es binaria, asignamos a una categoría el valor de 1 y a otra el valor de 0.
Asumiremos que la variable categórica es numérica para los efectos de todo cálculo.
Sin embargo, la interpretación se centrará en la comparación de categorías 0 y 1, nunca se interperará valores intermedios porquen o existen.
Agenda
Introducción al modelado de regresión
Modelo de Regresión Lineal Simple
Regresión Lineal Simple en R
Regresión Lineal Múltiple
Regresión Lineal Múltiple en R
Evaluación de Supuestos
El Modelo Lineal Generalizado
La regresión de Poisson
La regresión Binomial Negativa
Tablas de regresión reproducibles con {gtsummary}
Se usa la función lm() de R base. Sin embargo, la salida de esta no es muy informativa:
El modelo puede guardarse para realizar más operaciones sobre este. Por ejemplo, mejorar la salida.
Podemos usar summary() para ver resultados detallados.
Call:
lm(formula = y ~ x1, data = datos)
Residuals:
Min 1Q Median 3Q Max
-3.8666 -1.1168 -0.3487 1.3100 4.1336
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.06666 0.37316 -0.179 0.859
x1Tratamiento Nuevo 4.27094 0.52773 8.093 1.59e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.866 on 48 degrees of freedom
Multiple R-squared: 0.5771, Adjusted R-squared: 0.5683
F-statistic: 65.5 on 1 and 48 DF, p-value: 1.594e-10
Call:
lm(formula = y_peso_final ~ x3_peso_inicial, data = datos2)
Residuals:
Min 1Q Median 3Q Max
-10.0568 -4.7717 -0.8704 5.1824 10.4953
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.4317 2.6574 -2.044 0.0412 *
x3_peso_inicial 1.3447 0.1766 7.615 6.1e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.535 on 998 degrees of freedom
Multiple R-squared: 0.05491, Adjusted R-squared: 0.05397
F-statistic: 57.99 on 1 and 998 DF, p-value: 6.1e-14
\[y\_pesofinal = -5.4317 + 1.3447*x3\_pesoinicial + \epsilon_i\]
\[\epsilon_i \sim Normal(0, 5.535^2)\]
{broom} y su función tidy() podemos obtener también los intervalos de confianza:# A tibble: 2 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -5.43 2.66 -2.04 4.12e- 2 -10.6 -0.217
2 x3_peso_inicial 1.34 0.177 7.62 6.10e-14 0.998 1.69
Interpretación:
\(\beta_0\) o intercepto: Este viene a ser el valor promedio de \(y\) cuando todos los valores de \(x\) son 0. En este caso, cuando el peso inicial es cero kg. ¿Esto es posible?, por tal motivo, no se suele interpretar este valor.
\(\beta_1\) o coeficiente de regresión de x3_peso_inicial: Por cada 1 kg adicional de peso inicial, el valor promedio del peso final aumenta 1.43 kg (IC95% 1.00 a 1.69; p < 0.001).
Call:
lm(formula = y_peso_final ~ x1_tto, data = datos2)
Residuals:
Min 1Q Median 3Q Max
-7.7043 -1.6644 -0.0095 1.5849 8.5658
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.8771 0.1112 178.67 <2e-16 ***
x1_ttoTratamiento Nuevo -10.2325 0.1573 -65.04 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.488 on 998 degrees of freedom
Multiple R-squared: 0.8091, Adjusted R-squared: 0.8089
F-statistic: 4230 on 1 and 998 DF, p-value: < 2.2e-16
# A tibble: 2 × 7
term estimate std.error statistic p.value conf.low conf.h…¹
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 19.9 0.111 179. 0 19.7 20.1
2 x1_ttoTratamiento Nuevo -10.2 0.157 -65.0 0 -10.5 -9.92
# … with abbreviated variable name ¹conf.high
Interpretación:
\(\beta_0\) (Intercept): A menudo no se interpreta. Es el valor promedio de \(y_i\) cuando los valores de \(x\) son cero. En este caso, cuando el tratamien es cero (placebo). ¿Esto es posible?, sí es posible pero no es de ayuda para modelos explicativos, por lo que no se interpreta.
\(\beta1\) x1Tratamiento Nuevo: El promedio de peso final en quienes recibieron el tratamiento nuevo fue 10.23 kg menor que el de quienes recibieron placebo (Dif. medias = -10.23; IC95% -10.54 a -9.92; p < 0.001).
Agenda
Introducción al modelado de regresión
Modelo de Regresión Lineal Simple
Regresión Lineal Simple en R
Regresión Lineal Múltiple
Regresión Lineal Múltiple en R
Evaluación de Supuestos
El Modelo Lineal Generalizado
La regresión de Poisson
La regresión Binomial Negativa
Tablas de regresión reproducibles con {gtsummary}
El modelo de regresión lineal múltiple generaliza la RLS permitiendo evaluar la relación de varias covariables explicativas \(x\) sobre \(y_i\).
Componente sistemático:
\[E[Y | X_1 = x_{1i}, ..., X_p = x_{pi}] = E[y_i] = \mu_i = \beta_0 + \beta_1x_{1i} + ... + \beta_px_{pi}\]
Componente aleatoria:
\[y_i \sim N(\beta_0 + \beta_1x_{1i} + ... + \beta_px_{pi}, I\sigma^2)\]
\[\epsilon_i \sim N(0, \sigma^2) \]
Genera un hiperplano recto.
No podemos imaginarnos una imagen de esto, pero sí podemos analizarlo a nivel estadístico.
Agenda
Introducción al modelado de regresión
Modelo de Regresión Lineal Simple
Regresión Lineal Simple en R
Regresión Lineal Múltiple
Regresión Lineal Múltiple en R
Evaluación de Supuestos
El Modelo Lineal Generalizado
La regresión de Poisson
La regresión Binomial Negativa
Tablas de regresión reproducibles con {gtsummary}
Call:
lm(formula = y_peso_final ~ x1_tto, data = datos2)
Residuals:
Min 1Q Median 3Q Max
-7.7043 -1.6644 -0.0095 1.5849 8.5658
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.8771 0.1112 178.67 <2e-16 ***
x1_ttoTratamiento Nuevo -10.2325 0.1573 -65.04 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.488 on 998 degrees of freedom
Multiple R-squared: 0.8091, Adjusted R-squared: 0.8089
F-statistic: 4230 on 1 and 998 DF, p-value: < 2.2e-16
El modelo RLS solo incluye una covariable.
El modelo de RLM incluye 2 o más covariables.
Estas se agregan con un símbolo +
Notar que no se reportan intervalos de confianza al 95%.
Call:
lm(formula = y_peso_final ~ x1_tto + x3_peso_inicial, data = datos2)
Residuals:
Min 1Q Median 3Q Max
-5.5598 -1.4213 0.1343 1.0768 5.4482
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.94719 0.99689 -0.95 0.342
x1_ttoTratamiento Nuevo -10.25530 0.13111 -78.22 <2e-16 ***
x3_peso_inicial 1.38755 0.06614 20.98 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.073 on 997 degrees of freedom
Multiple R-squared: 0.8676, Adjusted R-squared: 0.8673
F-statistic: 3266 on 2 and 997 DF, p-value: < 2.2e-16
El paquete {broom} tiene funciones que facilitan obtener varios estadísticos de interés de los modelos de regresión.
La función tidy() permite obtener intervalos de confianza y otras medidas de interés.
El paquete {broom} tiene funciones que facilitan obtener varios estadísticos de interés de los modelos de regresión.
La función tidy() permite obtener intervalos de confianza y otras medidas de interés.
Además, retorna un tibble() que puede manipularse y luego embeberse en una tabla.
Primero se carga el paquete:
tibble() que puede ser manipulado, exportado a excel, y convertido a tabla.# A tibble: 3 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -0.947 0.997 -0.950 3.42e- 1
2 x1_ttoTratamiento Nuevo -10.3 0.131 -78.2 0
3 x3_peso_inicial 1.39 0.0661 21.0 3.10e-81
conf.int = TRUE mostramos también los intervalos de confianza.# A tibble: 3 × 7
term estimate std.error statistic p.value conf.low conf.…¹
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -0.947 0.997 -0.950 3.42e- 1 -2.90 1.01
2 x1_ttoTratamiento Nuevo -10.3 0.131 -78.2 0 -10.5 -10.0
3 x3_peso_inicial 1.39 0.0661 21.0 3.10e-81 1.26 1.52
# … with abbreviated variable name ¹conf.high
| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | -0.9471917 | 0.99689433 | -0.9501426 | 3.422701e-01 | -2.903444 | 1.009060 |
| x1_ttoTratamiento Nuevo | -10.2553009 | 0.13110928 | -78.2194924 | 0.000000e+00 | -10.512583 | -9.998019 |
| x3_peso_inicial | 1.3875541 | 0.06613671 | 20.9800905 | 3.097784e-81 | 1.257771 | 1.517337 |
mod %>%
tidy(conf.int = TRUE) %>%
mutate(estimate = round(estimate, 2),
conf.low = round(conf.low, 2),
conf.high = round(conf.high, 2),
p.value2 = case_when(
p.value < 0.001 ~ "<0.001",
p.value >= 0.001 ~ as.character(round(p.value, 3))
)) %>%
select(term, estimate, conf.low, conf.high, p.value2) %>%
gt()| term | estimate | conf.low | conf.high | p.value2 |
|---|---|---|---|---|
| (Intercept) | -0.95 | -2.90 | 1.01 | 0.342 |
| x1_ttoTratamiento Nuevo | -10.26 | -10.51 | -10.00 | <0.001 |
| x3_peso_inicial | 1.39 | 1.26 | 1.52 | <0.001 |
mod %>%
tidy(conf.int = TRUE) %>%
mutate(estimate = round(estimate, 2),
conf.low = round(conf.low, 2),
conf.high = round(conf.high, 2),
p.value2 = case_when(
p.value < 0.001 ~ "<0.001",
p.value >= 0.001 ~ paste("= ", round(p.value, 3))
)) %>%
mutate(
`Coeficiente (IC95%), p valor` =
glue("{estimate} (IC95% {conf.low} a {conf.high}), p {p.value2}"),
Variables = c("Intercepto", "Tratamiento nuevo vs. Placebo", "Peso inicial (kg)")
) %>%
select(Variables, `Coeficiente (IC95%), p valor`) %>%
gt()| Variables | Coeficiente (IC95%), p valor |
|---|---|
| Intercepto | -0.95 (IC95% -2.9 a 1.01), p = 0.342 |
| Tratamiento nuevo vs. Placebo | -10.26 (IC95% -10.51 a -10), p <0.001 |
| Peso inicial (kg) | 1.39 (IC95% 1.26 a 1.52), p <0.001 |
Call:
lm(formula = y_peso_final ~ x1_tto + x3_peso_inicial, data = datos2)
Residuals:
Min 1Q Median 3Q Max
-5.5598 -1.4213 0.1343 1.0768 5.4482
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.94719 0.99689 -0.95 0.342
x1_ttoTratamiento Nuevo -10.25530 0.13111 -78.22 <2e-16 ***
x3_peso_inicial 1.38755 0.06614 20.98 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.073 on 997 degrees of freedom
Multiple R-squared: 0.8676, Adjusted R-squared: 0.8673
F-statistic: 3266 on 2 and 997 DF, p-value: < 2.2e-16
\[y\_pesofinal = -0.94719 -10.25530*x1ttoTratamientoNuevo + 1.3875*x3\_pesoinicial + \epsilon_i\]
\[\epsilon_i \sim Normal(0, 2.073^2)\]
| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | -0.9471917 | 0.99689433 | -0.9501426 | 3.422701e-01 | -2.903444 | 1.009060 |
| x1_ttoTratamiento Nuevo | -10.2553009 | 0.13110928 | -78.2194924 | 0.000000e+00 | -10.512583 | -9.998019 |
| x3_peso_inicial | 1.3875541 | 0.06613671 | 20.9800905 | 3.097784e-81 | 1.257771 | 1.517337 |
\(\beta_0\) o intercepto: Este viene a ser el valor promedio de \(y\) cuando todos los valores de \(x\) son 0. En este caso, cuando el peso inicial es cero kg y cuando el tratamiento es placebo. ¿Esto es posible?, por tal motivo, no se suele interpretar este valor.
\(\beta_2\) o coeficiente de regresión de x1_ttoTratamiento Nuevo: El promedio de peso final en quienes recibieron el tratamiento nuevo fue 10.26 kg menor que el de quienes recibieron placebo, luego de ajustar por peso inicial (Dif. medias = -10.26; IC95% -10.51 a -9.99; p < 0.001).
\(\beta_1\) o coeficiente de regresión de x3_peso_inicial: Por cada 1 kg adicional de peso inicial, el valor promedio del peso final aumenta 1.39 kg, luego de ajustar por tatamiento recibido (IC95% 1.26 a 1.52; p < 0.001).
Agenda
Introducción al modelado de regresión
Modelo de Regresión Lineal Simple
Regresión Lineal Simple en R
Regresión Lineal Múltiple
Regresión Lineal Múltiple en R
Evaluación de Supuestos
El Modelo Lineal Generalizado
La regresión de Poisson
La regresión Binomial Negativa
Tablas de regresión reproducibles con {gtsummary}
Los errores (\(\epsilon_i\)) son medidas de la población a la que no tenemos acceso.
Los residuos (\(e_i\)) son el análogo a los errores pero obtenidos de la muestra observada.
Podemos usar los residuos para evaluar algunos supuestos sobre los errores.
Supuestos estadísticos del modelo
Linealidad
Independencia de observaciones
Homocedasticidad de los errores \(\epsilon_i\)
Normalidad de los errores \(\epsilon_i\) o de \(y_i\).
No problemas con la regresión:
Supuestos si queremos generalizar a una población finita bien definida
La muestra es representativa de la población.
Cuando no lo tenemos, solo podemos generalizar a una población que sabemos que existe pero no podemos definir. ¿Qué tan relevante puede ser esto?
Hay asignación aleatoria
Cuando no lo tenemos, tenemos que poder asumir (¿ingenuamente?) que se puede emular la asignación aleatoria de alguna manera:
En realidad, los supuestos de los modelos lineales son sobre el comportamiento probabilístico de \(y_i\).
Sin embargo, la idea de la existencia de los errores y de sus valores observados en la muestra, residuos resulta útil para evaluar supuestos.
Permiten reducir un problema de muchas dimensiones a solo 1 o 2 dimensiones.
Son como las placas radiográficas para el diagnóstico de los modelos.
Se usan los residuos para explorar el comportamiento de los \(y_i\) o los errores \(\epsilon\).
Preferiblemente usar gráficos de residuos.
Pruebas de hipótesis que usan residuos tienen los mismos problemas que discutimos en clases anteriores.
Podríamos usarlas para complementar análisis cuando los tamaños de muestra no son ni muy pequeños ni muy grandes.
La función check_model del paquete {performance} genera un panel de gráficos muy útil para evalur estos supuestos.
Podemos complentar el análisis de supuestos con funciones del paquete {car}.
Se puede evaluar si la homocedasticidad es consistente según cada variable predictora.
Si no lo es, se puede optar por modelar esta heterogeneidad de varianzas.
Se sugiere usar residuos estudentizados.
El supuesto de linealidad es sobre los coeficientes de regresión \(\beta\), no sobre las covariables.
Las variables X deben estar en una forma apropiada para que el supuesto se cumpla.
Es bien difícil que exista linealidad en la realida, pero puede ocurrir en raras y excepcionales ocasiones.
Se sugiere asumir no linealidad y pre-planear un modelamiento no lineal.
Entre los métodos que pueden usarse, tenemos:
Splines: Bastante usado y sugerido en bioestadística. Útil para ajustar por variables continuas.
Modelamiento Multivariablede polinomios fraccionales. También usado y recomendado en literatura biomédica. Útil para modelar forma como objetivo principal.
Polinomios. Menos flexible, puede ser útil si se conoce bien la relación o se busca mejorar ajuste.
Modelos aditivos generalizados. Útil si se buscar modelar la relación. Complejos y requieren muchos datos.
Veamos un ejemplo de modelamiento continuo con splines:
Evite categorizar la variable continua
Categorizar es muy malo: se pierde información y se corre el riesgo de sesgar resultados.
Si se quiere ajustar por variables continuas, use Splines o Polinomios fraccionales. No requiere interpretar sus resultados, pero si ajsutar bien!
Si se quiere evaluar la relación de la variable continua, planee un método estadístico para modelar la forma sin asumir linealidad.
Presuponga que la relación no es lineal.
Modelo y responda su pregunta. Si la relación es lineal, el modelo más complejo revelerá una línea recta.
No homogeneidad de varianzas
Podemos usar una estimación robusta de la varianza.
Los paquetes {sanwich} y {lmtest} proporcionan funciones útiles para esto.
Es bien difícil de creer que existe homogeneidad de varianzas en la vida real (salvo muy raras y excepcionales ocasiones).
# A tibble: 3 × 7
term estimate std.error statistic p.value conf.low conf.…¹
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -0.947 1.05 -0.906 3.65e- 1 -3.00 1.10
2 x1_ttoTratamiento Nuevo -10.3 0.131 -78.1 0 -10.5 -10.0
3 x3_peso_inicial 1.39 0.0692 20.0 2.45e-75 1.25 1.52
# … with abbreviated variable name ¹conf.high
Si distribución es normal (cosa que no podemos saber con certeza), podemos dejar de preocuparnos por este supuesto.
Si se cumple TLC, podemos dejar de preocuparnos por este supuesto.
Si no se cumple TLC o hay dudas razonables, podemos optar por alguna de las siguientes alternativas:
Agenda
Introducción al modelado de regresión
Modelo de Regresión Lineal Simple
Regresión Lineal Simple en R
Regresión Lineal Múltiple
Regresión Lineal Múltiple en R
Evaluación de Supuestos
El Modelo Lineal Generalizado
La regresión de Poisson
La regresión Binomial Negativa
Tablas de regresión reproducibles con {gtsummary}
Modelo lineal que permite modelar desenlaces de varios tipos de distribuciones.
Generaliza el modelo de regresión lineal.
Permite que \(Y_i\) siga otras distribuciones.
Componente sistemático:
\[ g(E(Y|x_{1i}, ..., x_{pi})) = g(E(Y_i)) = \eta_i = \beta_0 + \beta_1x_{1i} + ...+ \beta_px_{ip} \]
\(g()\) es la función de enlace.
\(\eta_i\) es el predictor linear.
\(E(Y|x_{1i}, ..., x_{pi}) = \mu_i\)
Componente aleatorio:
\[ Y_i \sim Distribucion~de~la~Familia~Exponencial \]
| Variable respuesta | Distribución de FE | Función de enlace canónica \(g()\) | Otras funciones de enlace comunes |
|---|---|---|---|
| Binaria | Bernoulli (Binomial con n = 1) | \(logit()\) | \(log()\) |
| Conteo | Binomial (con n > 1) | \(logit()\) | \(log()\) |
| Poisson | \(log()\) | ||
| Binomial negativo | \(log(\mu + k)\) | \(log()\) | |
| Continua positiva | Gamma | \(\frac{1}{\mu}\) | |
| Gausiana inversa |
Hace uso de estimación de máxima verosimilitud (MV).
Salvo el caso normal (donde MV = MCO), no existe solución cerrada para obtener los estimadores de MV.
métodos numéricos: Fisher Scoring, Newton Raphson, etc.No siempre la función de verosimilitud tiene un máximo.
Solo cuando se usa la función de enlace canónica.
Caso contrario, puede no tener solución única y hay problemas de convergencia.
Agenda
Introducción al modelado de regresión
Modelo de Regresión Lineal Simple
Regresión Lineal Simple en R
Regresión Lineal Múltiple
Regresión Lineal Múltiple en R
Evaluación de Supuestos
El Modelo Lineal Generalizado
La regresión de Poisson
La regresión Binomial Negativa
Tablas de regresión reproducibles con {gtsummary}
\[y_i \sim Poisson(\beta_0 + \beta_1x_{1i} + ...+ \beta_px_{ip})\]
\[log(E(y_i)) = \eta_i\]
\[\eta_i = \beta_0 + \beta_1x_{1i} + ...+ \beta_px_{ip}\]
\[y_i \sim Poisson(\eta_i)\]
\(log()\) es la función de enlace canónica: solución única para MV y no problemas de convergencia por esto.
Si usamos la función identidad de la regresión lineal, el modelo quedaría planetado de esta manera:
\[E(y_i) = \beta_0 + \beta_1x_{1i} + ...+ \beta_px_{ip}\]
razón de medias, medida interpretable.Porque la distribución de \(y_i\) no es normal, es una variable de conteo.
El principal problema de esto, es que al ser Poisson, la \(media = varianza = \lambda\), por lo que a mayor valor de la media, la varianza aumentará.
Lo que implica que \(y_i\) es una v.a. heterocedástica.
El modelo normal necesita homocedasticidad, caso contrario, tiene que corregirse de alguna manera.
Poisson no necesita esto, su modelo es heterocedastico por naturaleza, lo que hace más eficiente la estimación:
La regresión de Poisson permite retornar directamente razón de medias (RM).
Los coeficientes de regresión \(\beta\) del modelo son \(log(RM)\), por lo tanto, podemos exponenciarlos para obtener los OR:
\[\beta = log(RM)\]
\[e^\beta = RM\]
Call:
glm(formula = docvis ~ female + age, data = md_visit)
Deviance Residuals:
Min 1Q Median 3Q Max
-5.139 -2.792 -1.623 0.593 118.711
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.442309 0.168002 -2.633 0.00848 **
female 0.980325 0.082751 11.847 < 2e-16 ***
age 0.071883 0.003682 19.522 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 33.29453)
Null deviance: 671518 on 19608 degrees of freedom
Residual deviance: 652772 on 19606 degrees of freedom
AIC: 124390
Number of Fisher Scoring iterations: 2
Se especifica la ecuación.
Call:
glm(formula = docvis ~ female + age, family = poisson, data = md_visit)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.3562 -2.2285 -1.1346 0.3396 26.8860
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.0368083 0.0176597 -2.084 0.0371 *
female 0.3092910 0.0081168 38.105 <2e-16 ***
age 0.0227500 0.0003624 62.768 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 122270 on 19608 degrees of freedom
Residual deviance: 116363 on 19606 degrees of freedom
AIC: 153636
Number of Fisher Scoring iterations: 6
Se indica la familia de distribución de \(y_i\).
Call:
glm(formula = docvis ~ female + age, family = poisson(link = "log"),
data = md_visit)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.3562 -2.2285 -1.1346 0.3396 26.8860
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.0368083 0.0176597 -2.084 0.0371 *
female 0.3092910 0.0081168 38.105 <2e-16 ***
age 0.0227500 0.0003624 62.768 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 122270 on 19608 degrees of freedom
Residual deviance: 116363 on 19606 degrees of freedom
AIC: 153636
Number of Fisher Scoring iterations: 6
Se indica la familia de distribución de \(y_i\).
mod <- glm(docvis ~ female + age,
family = poisson(link = "identity"),
data = md_visit)
summary(mod)
Call:
glm(formula = docvis ~ female + age, family = poisson(link = "identity"),
data = md_visit)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.1714 -2.2486 -1.1320 0.3312 26.8394
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.18852 0.04629 -4.073 4.64e-05 ***
female 1.00553 0.02510 40.056 < 2e-16 ***
age 0.06581 0.00110 59.848 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 122270 on 19608 degrees of freedom
Residual deviance: 116471 on 19606 degrees of freedom
AIC: 153743
Number of Fisher Scoring iterations: 6
mod <- glm(docvis ~ female + age,
family = poisson(link = "log"),
data = md_visit)
mod %>%
tidy() %>%
gt()| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -0.03680829 | 0.0176596601 | -2.084315 | 0.03713153 |
| female | 0.30929097 | 0.0081167612 | 38.105220 | 0.00000000 |
| age | 0.02275001 | 0.0003624464 | 62.767938 | 0.00000000 |
log(razón de medias), para volverlos interpretables, hay que aplicar antilogaritmo: exp().| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 0.9638609 | 0.0176596601 | -2.084315 | 0.03713153 |
| female | 1.3624587 | 0.0081167612 | 38.105220 | 0.00000000 |
| age | 1.0230108 | 0.0003624464 | 62.767938 | 0.00000000 |
exponentiate = TRUE.| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | 0.9638609 | 0.0176596601 | -2.084315 | 0.03713153 | 0.9310334 | 0.9977672 |
| female | 1.3624587 | 0.0081167612 | 38.105220 | 0.00000000 | 1.3409630 | 1.3843148 |
| age | 1.0230108 | 0.0003624464 | 62.767938 | 0.00000000 | 1.0222845 | 1.0237380 |
Call:
glm(formula = docvis ~ female + age, family = poisson(link = "log"),
data = md_visit)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.3562 -2.2285 -1.1346 0.3396 26.8860
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.0368083 0.0176597 -2.084 0.0371 *
female 0.3092910 0.0081168 38.105 <2e-16 ***
age 0.0227500 0.0003624 62.768 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 122270 on 19608 degrees of freedom
Residual deviance: 116363 on 19606 degrees of freedom
AIC: 153636
Number of Fisher Scoring iterations: 6
female: El número medio de visitas anuales al médico en mujeres fue 20% veces más el de los hombres (RM = 1.33; IC95% 1.31 a 1.35; p < 0.001)
age: Por cada incremento de la edad en un año, el número medio de visitas anuales al médico se incrementa en 1% (RM = 1.017; IC95% 1.016 a 1.018; p < 0.001).
Linealidad del \(log(y_i)\) respecto a la combinación lineal de predictores.
Observaciones son independientes.
\(Y_i\) sigue distribución de Poisson.
No problemas de regresión:
No puntos influyentes
No colinealidad: Solo cuando esta es un problema.
Supuestos específicos si se busca generalizar a poblaciones conocidas, hacer inferencias causales o ambas.
Agenda
Introducción al modelado de regresión
Modelo de Regresión Lineal Simple
Regresión Lineal Simple en R
Regresión Lineal Múltiple
Regresión Lineal Múltiple en R
Evaluación de Supuestos
El Modelo Lineal Generalizado
La regresión de Poisson
La regresión Binomial Negativa
Tablas de regresión reproducibles con {gtsummary}
\(y_i|\lambda_i \sim Poisson(\lambda_i)\) y \(\lambda_i \sim Gamma(\mu_i, \psi)\)
Entonces, es posible mostrar que \(y_i\) sigue una distribución binomial negativa.
Asimismo, el modelo:
\[y_i \sim BN(\beta_0 + \beta_1x_{1i} + ...+ \beta_px_{ip}, \psi)\]
\[log(E(y_i)) = \eta_i\]
\[\eta_i = \beta_0 + \beta_1x_{1i} + ...+ \beta_px_{ip}\]
\[y_i \sim BN(\eta_i, \psi)\]
Porque no siempre la variable \(y_i\) seguirá una distribución de Poisson.
Si sigue una distribución BN, entonces la varianza es mayor a la media (sobredispersion).
La estimación del error estándar deberá tener esto en cuenta.
Caso contrario, sería inválido en dirección anticorservadora.
La función de enlace de la BN no es log(), tampoco identiy().
Sin embargo, se prefiere usar log() para obtener resultados interpretables: razón de medias.
Esto conlleva un problema, no siempre hay convergencia:
- Cuando hay sobredispersión, regresión binomial negativa no siempre es la opción factible.
- Otra opción puede ser usar regresión quasipoisson o un estimador de varianza robusta.La regresión BN permite retornar directamente razón de medias (RM).
Los coeficientes de regresión \(\beta\) del modelo son \(log(RM)\), por lo tanto, podemos exponenciarlos para obtener las RM:
\[\beta = log(RM)\]
\[e^\beta = RM\]
Call:
glm.nb(formula = docvis ~ female + age, data = md_visit, init.theta = 0.4838144807,
link = log)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.5688 -1.3167 -0.4629 0.1127 6.5046
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.0366979 0.0457077 -0.803 0.422
female 0.3370406 0.0222361 15.157 <2e-16 ***
age 0.0224236 0.0009913 22.621 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(0.4838) family taken to be 1)
Null deviance: 21019 on 19608 degrees of freedom
Residual deviance: 20224 on 19606 degrees of freedom
AIC: 85857
Number of Fisher Scoring iterations: 1
Theta: 0.48381
Std. Err.: 0.00659
2 x log-likelihood: -85849.42400
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -0.03669795 | 0.0457077355 | -0.8028827 | 4.220425e-01 |
| female | 0.33704062 | 0.0222361245 | 15.1573453 | 6.775303e-52 |
| age | 0.02242362 | 0.0009912903 | 22.6206366 | 2.715341e-113 |
log(razón de medias), para volverlos interpretables, hay que aplicar antilogaritmo: exp().| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 0.9639673 | 0.0457077355 | -0.8028827 | 4.220425e-01 |
| female | 1.4007960 | 0.0222361245 | 15.1573453 | 6.775303e-52 |
| age | 1.0226769 | 0.0009912903 | 22.6206366 | 2.715341e-113 |
exponentiate = TRUE.| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | 0.9639673 | 0.0457077355 | -0.8028827 | 4.220425e-01 | 0.8824872 | 1.053230 |
| female | 1.4007960 | 0.0222361245 | 15.1573453 | 6.775303e-52 | 1.3412238 | 1.463041 |
| age | 1.0226769 | 0.0009912903 | 22.6206366 | 2.715341e-113 | 1.0207472 | 1.024612 |
Call:
glm.nb(formula = docvis ~ female + age, data = md_visit, init.theta = 0.4838144807,
link = log)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.5688 -1.3167 -0.4629 0.1127 6.5046
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.0366979 0.0457077 -0.803 0.422
female 0.3370406 0.0222361 15.157 <2e-16 ***
age 0.0224236 0.0009913 22.621 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(0.4838) family taken to be 1)
Null deviance: 21019 on 19608 degrees of freedom
Residual deviance: 20224 on 19606 degrees of freedom
AIC: 85857
Number of Fisher Scoring iterations: 1
Theta: 0.48381
Std. Err.: 0.00659
2 x log-likelihood: -85849.42400
# A tibble: 3 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0.964 0.0457 -0.803 4.22e- 1 0.882 1.05
2 female 1.40 0.0222 15.2 6.78e- 52 1.34 1.46
3 age 1.02 0.000991 22.6 2.72e-113 1.02 1.02
female: El número medio de visitas anuales al médico en mujeres fue 20% veces más el de los hombres (RM = 1.40; IC95% 1.34 a 1.46; p < 0.001)
age: Por cada incremento de la edad en un año, el número medio de visitas anuales al médico se incrementa en 2.3% (RM = 1.017; IC95% 1.021 a 1.025; p < 0.001).
Linealidad del \(log(y_i)\) respecto a la combinación lineal de predictores.
Observaciones son independientes.
\(Y_i\) sigue distribución de Poisson.
No problemas de regresión:
No puntos influyentes
No colinealidad: Solo cuando esta es un problema.
Supuestos específicos si se busca generalizar a poblaciones conocidas, hacer inferencias causales o ambas.
Agenda
Introducción al modelado de regresión
Modelo de Regresión Lineal Simple
Regresión Lineal Simple en R
Regresión Lineal Múltiple
Regresión Lineal Múltiple en R
Evaluación de Supuestos
El Modelo Lineal Generalizado
La regresión de Poisson
La regresión Binomial Negativa
Tablas de regresión reproducibles con {gtsummary}
Podemos usar la librería {gtsummary} para esto.
Veamos un ejemplo.
Podemos reportar la tabla de regreion multivarible de la siguiente manera:
# A tibble: 3 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 11.0 0.0505 218. 0 10.9 11.1
2 age 0.0110 0.000341 32.3 5.81e-218 0.0104 0.0117
3 sex -0.474 0.0258 -18.3 8.44e- 74 -0.524 -0.423
tbl_regression() de {gtsummary}:tbl_uvregression() de {gtsummary}:tbl_merge():tabla_final <- tbl_merge(list(tabla_univ, tabla_multi), tab_spanner = c("Modelos crudos", "Modelo ajustado"))
tabla_final| Characteristic | Modelos crudos | Modelo ajustado | |||||
|---|---|---|---|---|---|---|---|
| N | Beta | 95% CI1 | p-value | Beta | 95% CI1 | p-value | |
| age | 10,000 | 0.01 | 0.01, 0.01 | <0.001 | 0.01 | 0.01, 0.01 | <0.001 |
| sex | 10,000 | -0.47 | -0.53, -0.42 | <0.001 | -0.47 | -0.52, -0.42 | <0.001 |
| 1 CI = Confidence Interval | |||||||
https://github.com/psotob91
percys1991@gmail.com
R Aplicado a los Proyectos de Investigación - Sesión 8